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Abstract 

We carry out a finite density calculation based on a canonical approach which is designed to address the overlap problem. 
Two degenerate flavor simulations are performed using Wilson gauge action and Wilson fermions on 4'* lattices, at temperatures 
close to the critical temperature Tc « 170 MeV and large densities (5 to 20 times nuclear matter density). In this region, we 
find that the algorithm works well. We compare our results with those from other approaches. 
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I. INTRODUCTION 

The phase structure of QCD at finite temperature and 
finite density is relevant for a variety of phenomena: from 
subtle modifications of cross-section in high energy colli- 
sions of nuclei, to exotic states of nuclear matter in neu- 
tron stars. Due to asymptotic freedom we can use pertur- 
bation theory to study the quark-gluon plasma at suffi- 
ciently large temperatures. However, the regions of inter- 
est for heavy ion collision experiments and astrophysics 
are essentially non-perturbative. Numerical studies of 
QCD are extremely helpful in providing a quantitative 
understanding of the phase structure in these regions. 

At zero baryon density, it has been known for quite 
some time that QCD undergoes a transition from a 
confined phase to a deconfined phase at a temperature 
Tc « 170 MeV. Lattice QCD suggests that the transition 
is in fact a smooth crossover. This is expected to turn 
into a first order phase transition as the baryon density 
is increased. A schematic picture of the expected phase 
diagram is presented in Figure ^ 

The position of the second order transition point, 
where the crossover turns into a first order phase transi- 



quark-gluon plasma 



hadronic phase 



FIG. 1: Schematic picture of the expected QCD phase dia- 
gram. The solid line represents a first order phase transition, 
the dot a second order phase transition and the dashed line 
represents a crossover. 



tion, is very important in providing a quantitative under- 
standing of the QCD phase diagram. Close to the critical 
temperature the relevant degrees of freedom are the glu- 
ons and three flavors of quarks, the light quarks, up and 
down, and the strange quark. The shape of the curve 
seems to depend very little on the masses of the quarks, 
but the position of the second order transition point de- 
pends strongly on the mass of the strange quark. All nu- 
merical simulations treat the light quarks as degenerate. 
If the strange quark mass is taken to be equal to the mass 
of the light quarks, we have a theory with three degen- 
erate flavors. In this case, for low enough quark masses, 
the zero density phase transition is expected to be first 
order, the second order point disappears. As the strange 
quark mass is increased, the zero density phase transi- 
tion becomes a crossover, and the second order phase 
transition point moves to larger and larger densities. As 
the strange quark becomes infinitely heavy, only the light 
quarks remain dynamically relevant. The position of the 
second order phase transition point is not the same as in 
the physically relevant case; however, qualitatively the 
picture remains the same. This is why two degenerate 
flavor QCD is interesting as a testbed for methods to 
simulate finite density QCD. 

Simulations at finite temperature and zero baryon den- 
sity can be performed using standard lattice techniques. 
However, non-zero baryon density calculations remain 
one of the challenges of Lattice QCD. The reason is that, 
at non-zero chemical potential, the fermionic determi- 
nant becomes complex and the standard Monte Carlo 
methods fail since the integrand is no longer real and 
positive definite. The usual approach is to split the inte- 
grand in two parts, one that is real and positive and can 
be employed to generate an ensemble of configurations, 
and another one that includes the complex phase of the 
determinant and is folded into the observables. For clar- 
ity, let's write the grand canonical partition function for 
Lattice QCD: 



-S,{U)-Si(p.-U,-4>,^) 



(1) 



The fermionic part of the action, 

Sf{m,fj,;U,'ip,ip) = 'tpM{m,^i;U)ip, (2) 
is a quark bilinear and we can perform the path integra- 
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tion analytically, 



Z{V,T,fi) - y"l?Ue~^»(^)[|detM(m„M,;C/), 



(3) 



where M is the quark matrix and and /i^ are the mass 
and the chemical potential for flavor i. The gluonic part 
of the integrand, e"'^"'-^-', is real and positive, whereas 
the fermionic part is only guaranteed to be real when the 
chemical potential is zero. In the case of two degenerate 
flavors, after setting Hi — ^2 — IJ^, the partition function 
becomes 

Z{V,T,fi)^ / PUe-^«(^' detM(TO,/i;[/)2. (4) 



The standard approach, the Glasgow reweighting method 
01 , is to split the fermionic part into a real positive part 
and a phase factor 



det M{fif ^ det M{fi = 0)^ x 



detAf(^)^ 
detM(^ = 0)2^ 



(5) 



where we dropped some redundant indices. We can then 
apply the standard Monte Carlo techniques to generate 
an ensemble according to the measure. 



P{U) = e"^s(c^) det M(/i = 0)^ 
and then insert the phase factor. 



0{fi;U) = 



into the observable 



det M{p;U)^ 
det M{p = 0;[/)2' 



{0{U)e{^l;U))^^, 



(6) 



(7) 



(8) 



To establish some terminology, we will call the ensemble 
generated with the weight P{U) the generated ensemble 
and, in a manner of speaking, we will be calling target en- 
semble the one that would be generated using the weight 
derived from the true action. For the second term, the 
word ensemble is used loosely since for a complex inte- 
grand the concept of ensemble is, at best, ambiguous. 

There are two major problems with the reweighting 
approach: the sign problem and the overlap problem. 
The sign problem appears when the phase factor, 9{fj.; U), 
averages to a value too close to zero on the generated 
ensemble. By close, we mean an average value that is 
smaller than the error. In that case, all the measurements 
will have sizeable error bars and the method fails since 
we need extremely large ensembles to get reasonable error 
bars. 

The second problem appears when the generated en- 
semble and the target ensemble overlap poorly; for ex- 
ample, they might be in different phases. More precisely, 
take an observable (in our example, the order parameter 
for the phase transition) , if it happens that the histogram 



of this observable in the generated ensemble overlaps very 
poorly with the histogram in the target ensemble, then 
the value that we get via reweighting will be wrong. This 
problem is more serious than the sign problem since there 
is no indication when the measurement fails; the error 
bars can be deceptively small P, Q . 



Recently, a lot of progress has been made in studying 
the phase diagram at temperatures around Tc and small 
chemical potential 0, 0, B IE 0- The reason is that, 
in this region, the sign problem is manageable and the 
overlap problem is expected to be under control. The 
methods employ a more or less sophisticated form of 
reweighting ^ |E IS or some form of analytical contin- 
uation from imaginary chemical potential J^. The main 
results are the shape of the phase transition curve around 
Tc and the location of the second order phase transition 
point for quark masses close to the physical masses. All 
these simulations seem to be free of the sign problem, 
but it is not clear whether the overlap problem is indeed 
under control. One way to make sure that the results 
are correct is to either use methods that are proved to 
be free of overlap problem, or methods that are different 
enough but produce the same results. For small values 
of ^/Tc <C 1, different methods seem to agree. However, 
there is only one result for the location of the second or- 
der phase transition point, and it occurs at rather large 
value of the chemical potential. It is thus important to 
ask whether this result is reliable. 



In light of the problems mentioned above, it is imper- 
ative that new methods be developed to simulate QCD 
at finite density. All the methods mentioned above are 
based on the grand canonical partition function. Far 
fewer attempts have been made to simulate QCD using 
the canonical partition function [^I^ITo|. In this paper, 
we will present simulations based on a method that em- 
ploys the canonical partition function [ill IT3. IT^ . The 
main idea is that, to avoid the overlap problem, it is es- 
sential to generate an ensamble that is based on the pro- 
jected determinant, instead of reweighting. Moreover, to 
reduce the determinant fluctuations, the updating pro- 
cess is broken into two steps: an HMC proposal and an 
accept/reject step based on determinant ratios. These 
runs are exploratory in nature, using very small lattices 
and rather large quark masses. The main goal is to de- 
termine the feasibility of the algorithm and explore the 
available phase space. 



The paper is organized as follows: in section II we will 
introduce the canonical ensemble for QCD, in section III 
we present the algorithm we employed, in section IV we 
discuss the performance of the algorithm and in section 
V we present the physical results. We then conclude 
by attempting a physical interpretation of our results in 
section VI. 
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II. CANONICAL PARTITION FUNCTION 

The simplest way to show how to build the canonical 
ensemble in Lattice QCD is to start from the fugacity 
expansion, 



Z{V,T,t,) = Y,ZciV,T,n)e'^-/^, 

n 



(9) 



where n is the net number of quarks (number of quarks 
minus the number of anti-quarks) and Zc is the canonical 
partition function. We note here that on a finite lattice, 
the maximum net number of quarks is limited by the 
Pauli exclusion principle. Using the fugacity expansion, 
it is easy to see that we can write the canonical partition 
function as a Fourier transform of the grand canonical 
partition function, 

Zc{V,T,n) = — d0e-'"^Z(y,T,Ai)|;.=.0T. (10) 

We will now specialize to the case of two degenerate 
flavors. We use the grand canonical partition function in 
Eq. Q to get 

Zc{V,T,n) = / PUe"^«(^)det„M2(J7), (11) 



where we define 
det„Af2([/) 



2tt jo 



-incp 



det M(m,^; C/)^|^=j0t- 
(12) 

It is worth pointing out that the canonical partition func- 
tion defined above sums over configurations where the 
total net number of quarks, n = m + n2, is fixed. If we 
want to fix the net quark number for each flavor then we 
would use detmM{U) x dein^M{U). 

For our study, we will be using Wilson fermions. To 
introduce a non-zero chemical potential, the fcrmion ma- 
trix at zero chemical potential. 



[M{U)]x,y = 5x,y - At y^(l ~ -f^)U^{x)5x+[i,v 
4 

- ^^(l+7p)f/^(2/)<5.,,+^, (13) 

is altered 0^ by introducing a bias for time forward prop- 
agation in the hopping matrix. More specifically, the 
hopping in the time direction is altered 

(1 + 74)C/1(2/) ^ (l + 74)^7l(y)e^^ 
{l-li)Ui{x) ^ (1 - 74)^74(x)e-''^ (14) 



We can perform a change of variables 



^(f, Xi) 



^^{X, Xi) 



to restore the original form of the hopping matrix except 
on the last time slice. In terms of these new variables, we 
can write the fermionic matrix as M(Utf,) — M{m,fj.; U) 
where M is defined in Eq. (|13|l and 



X4, = Nt,iy = A 
otherwise. 



(16) 



This should not be viewed as a change of the gauge field 
variables but rather as a convenient way to write the 
fermionic matrix. 

In order to evaluate numerically the partition function 
in Eq. ((TT|l . we need to replace the continuous Fourier 
transform in Eq. H12|l with a discrete one. We will then 
redefine the projected determinant. 



N-l 



det„A/2({/) = ^J2 ^"""^^ det M{U^^y 



(17) 



where cj)j = and the parameter N defines the discrete 
Fourier transform. In the limit — > 00, we recover the 
original projected determinant. For finite N the partition 
function 

Zc{V,T,n) = [ VVe-^^^^'>(ktnM^{U), (18) 



will only be an approximation of the canonical partition 
function. Using the fugacity expansion we can show that 

00 

Zc{V,T,n)^ Zc{V,T,n + mN). (19) 

m— — 00 

If N and n are chosen such that \n + mN\ is minimal for 
TO = then Zc should be a good approximation to Zc 
as long as 



Zc{V,T,n + mN) 
Zc{V,T,n) 



«1, 



(20) 



for all m ^ 0. To understand better this condition take 
< n < N/2\ the largest contamination comes from 
Zc{V, T,N ~n). The ratio above is 



Zc-(u,r,iv-7 
Zc(y,r,n) 



F {V ,T ,N - n) - F {V ,T ,n) 



(21) 



ijj'{x,Xi) = e^''^Xf,a;4) 



(15) 



where F is the free energy of the system. For low tem- 
peratures, we expect that F(y,T, n) oc e~*^«l"l/-^, where 
Mb is the mass of the baryon. We see then that the 
approximation will hold as long as the temperature is 
low enough or the baryon mass is high. This assumption 
needs to be checked in our simulations; if it fails then we 
need to increase N . 



III. ALGORITHM 

In this section, we will present the algorithm we employ 
to simulate the partition function Zc- Directly simulat- 
ing the projected determinant in Eq. H17|) is known to 
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Projected determinant fluctuations Determinant ratio fluctuations 




50 100 150 200 250 50 100 150 200 250 



configuration configuration 

FIG. 2: Fluctuations of the fermionic part of the measure and of the accept /reject factor i^{U), defined in Eq. 13011 . as measured 
on an ensemble generated at /3 = 5.2 and n = 3. In both figures, we subtracted the average value so that the plots are centered 
around zero. 



face a fluctuation problem ^T3\ , since In det M = Tr In Af 
and TrlnM is proportional to the lattice volume. To 
alleviate the problem we split the Markov process in 
two steps: a proposal step based on HMC and an ac- 
cept/reject step based on the ratio of the projected deter- 
minant to the determinant used in the HMC step. Since 
the accept/reject is based on the determinant ratio, the 
fluctuations should be reduced and the acceptance rate 
enhanced. We shall test this numerically. 



A. Target measure 

To evaluate Zc using Monte Carlo techniques, we need 
the integrand to be real and positive. Using 75 hermitic- 
ity of the fermionic matrix, i.e. 



(22) 



we can easily prove that det M{U^) is real. This implies 
that 



detnM^(U) ^ (det_„M2 ([/))' 



(23) 



but it doesn't imply that the projected determinant 
is real. det„M^(t/) is real only if det M{U^) — 
detM([/_0), which is not true configuration by config- 
uration. We can prove that for a charge conjugation 
symmetric action, this property is true when averaged 
over the ensemble, i.e. [det M{IJ^)) = {det M{U-^)). In 
fact, using charge conjugation symmetry of the action, 
we can prove that 

Zc{V,T,n)^ Zc{V,T,-n). (24) 

This property allows us to rewrite the partition function 

Zc{V,T,n) = /pUe"^«(^)Redet„M2(f/). (25) 



Now the integrand is real but not necessarily positive. 
For the sake of the argument, let's set aside for a while 



the fact that the integrand may be negative and assume 
that the above expression can be evaluated using stan- 
dard Monte Carlo techniques. Even then, the fact that 
we can write the partition function using a real integrand 
is not sufficient. The goal of any simulation is to compute 
different observables, the partition function itself is not 
of much interest. If we are only interested in observables 
that are even under charge conjugation, then we could 
use the ensemble generated by the above action to com- 
pute them. For observables that are odd under charge 
conjugation an additional step is necessary: we have to 
reintroduce a phase. 

We want to emphasize that, if the above integrand is 
positive, the observables which are even under charge 
conjugation could be evaluated directly on the ensemble 
generated by the above action. Thus, we would have 
no reweighting involved, and no overlap problem. The 
observables that are odd under charge conjugation are 
not guaranteed to behave as well, but we assume that 
their behavior would be similar. At the worst, the extra 
phase might introduce a sign problem. 

We come back now to the positivity question. In the 
case that the integrand is not positive, we are forced to 
use the absolute value of the integrand as measure for our 
generated ensemble. The algorithm will then be designed 
to generate an ensemble according to the weight 



W{U) 



-Sc,{U} 



Redet„M2([/) 



(26) 



The sign will be folded into the observables. For a generic 
observable the sign will turn out to be some complex 
phase. 



aiU) 



det„M2(t/) 



Redet„M2(t/) 



(27) 



but for observables even under charge conjugation it will 
be just the sign of Redet„M^(J7). 
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From the above discussion, it is clear that as long as 
we don't have a sign problem, the results of our simula- 
tion for observables invariant under charge conjugation 
arc reliable. For the other observables, the sign problem 
might be more severe, but we expect that they will not 
have an overlap problem. 



B. HMC update 

Turning to the algorithmic issues, our approach to gen- 
erating an ensemble with weight W{U) is to employ a 
Metropolis accept/reject method. In short, the method 
employs a generating mechanism that proposes new con- 
figurations with weight W'{U) and then an accept/reject 
step is used to correct for the target weight. Ideally, the 
proposal mechanism would propose configurations with 
the weight W{U); in that case all new proposals will be 
accepted. In practice, it isn't always possible to design 
efficient proposal mechanism for every weight. The gen- 
eral approach is to use an efficient proposal mechanism 
to generate a weight W'{U) close to the target weight 
W{U). If successful, the acceptance rate would be high 
and the algorithm would be efficient. 

One possible solution is to use a heat-bath method to 
propose new configurations based on the weight W'{U) = 
g-5'g(c/)_ However, such an updating strategy would be 
inefficient since the fermionic part is completely disre- 
garded in the proposal step. The determinant, being an 
extensive quantity, can fluctuate wildly from one config- 
uration to the next in the pure gauge updating process 
To reduce the fluctuations, it was suggested 
[l^ that we should employ an HMC algorithm for the 
proposal step. In this case 



W'{U) = e-^«(^) det M{uy' 



(28) 



We will then accept the new configuration U' with the 
probability 



Pacc = min{I,c^((7')/^^(f/)}, 
where u; is the ratio of the weights 

Reckt„M2([/) 



W{U) 
W'{U) 



detA/2([/) 



(29) 



(30) 



We expect that this proposal mechanism will be more 
efficient. Although the fermionic part of the measure 



Redet„M2([/) 



varies significantly from one configura- 



tion to the next, the determinant ratio 



ujiU) 



1 

N 



N-l 

cos{if)in)e 

3=0 



Tr{logM{U^.)-logM{U)) 



(31) 



is expected to fluctuate less. We base our expectation 
on the fact that, in the ratio, the leading fluctuations 









J^ — 
j 




1 







^ 



50 100 150 200 250 300 
iteration 

FIG. 3: Polyakov loop argument as a function of the simula- 
tion time. Note that toward the end the value is unchanged 
for almost 50 iterations. 



are removed by the Tr log difference of the quark matri- 
ces M{Uip ) and M{U). To check this, we compare the 
fluctuations in the fermionic part of the determinant and 
the determinant ratios in one of the ensembles generated 
in our runs. This is shown in Fig. |21 We see that the 
fluctuations are significantly reduced which results in a 
large boost in acceptance rate. 



C. Triality 

The canonical partition function, Eq. (|If |l . has a 
symmetry that is a direct consequence of the sym- 
metry of the grand canonical partition function at imag- 
inary chemical potential |l8|| . Under a transformation 
U ^ Utfi with (f) — ±27r/3, the gauge part of the action is 
invariant and 

det„Af2(C/) ^ det„M2([/±2./3) = e±''^"det„A/2({/). 

(32) 

We see then that when n is a multiple of 3 this transfor- 
mation leaves det„M^(t/) invariant. Consequently, the 
canonical partition action is invariant under this trans- 
formation. Incidentally, this symmetry of the gauge part 
of the action together with the transformation above of 
the fermionic part guarantees that the canonical parti- 
tion function will vanish when n is not a multiple of 3. 
However, this is no longer true if this symmetry is spon- 
taneously broken, which is the case in the deconfined 
phase. In this phase, there is no reason to expect that 
the canonical partition function should vanish when n is 
not a multiple of 3. 

The transformation rule above is preserved for the dis- 
crete case if we choose TV, the parameter that defines the 
Fourier transform, to be a multiple of 3. In our simu- 
lations, we will always choose N to satisfy this condi- 
tion. In this case, the remarks we made about Zc are 
valid for Zc and the projected determinant, det„Af^(L/) 
is invariant. Thus, the measure is symmetric under this 
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transformation, i.e. 

W{U) = W{U±2./z)- (33) 

However, the HMC weight, W'{U), does not have this 
symmetry since detM([/) is not invariant under this 
transformation. Because of this, our algorithm can be- 
come frozen for long periods of time. For example, in 
Fig. 121 we show how the argument of the Polyakov loop 
changes with the simulation time if we use the method 
presented so far. We notice that at the end of the simula- 
tion, when we tunnel to the sector where arg[P] « 27r/3, 
the update is frozen; the new proposals are rejected for 
a long time. This is due to the fact that HMC strongly 
prefers the sector. To understand this better, assume 
that we have a configuration Uq in the sector, where 
arg[P(J7o)] ~ 0, and denote with C7+ the configuration 
(t^o)27i-/3 with arg[P([/+)] w 27r/3. Then, we expect that 
det M2(C/o) > det NP{U+) since HMC prefers the sec- 
tor, but det„M^(?7o) = det„M^(J7+) since the projected 
determinant is symmetric under the transformations. 
Assume now that HMC proposes the accept/reject 
step will accept this since 

^{U+) ^ detM^([/o) 

Lu{Uo) detM2({/+) ^ • ^ ' 

However, in the next step HMC is likely to propose a new 
configuration in the sector since it favors it strongly. By 
the reverse of the argument above we have that 

u{U+) dctM2([/o) ^ ' 

and the new configuration will be very likely rejected. 

This means that although the algorithm will end up 
sampling the three sectors equally, as required by the 
symmetric weight W{U), two of the sectors will take a 
very long time to sample properly. To address this prob- 
lem, we introduce a hopping 0. Since the weight 
W(U) is symmetric under the Z^ transformation, we can 
intermix the regular updates with a change in the field 
variables U U±2tt/3- We will choose the sign randomly, 
with equal probability for each sign, to satisfy detailed 
balance. The new algorithm will sample all sectors in the 
same manner. 



IV. SIMULATION DETAILS 

Most of the computer time in these simulations is spent 
computing the determinant. There is a proposal that 
would employ a determinant estimator 19], but in this 
work we compute the determinant exactly using LU de- 
composition. This is a very expensive calculation con- 
sidering that even for the small lattices we used in this 
study the fermionic matrix has 3072 rows. Furthermore, 
the algorithm scales with the third power of the lattice 
four volume and it is not easily parallelizable. The high 



TABLE I: Simulation parameters. 





a(fm) 


m^(MeV) 


V ^(fm ^) 


T(MeV) 


5.00 


0.343(2) 


926(7) 


0.387(7) 


144(1) 


5.10 


0.322(4) 


945(13) 


0.468(17) 


153(2) 


5.15 


0.313(3) 


942(11) 


0.510(15) 


157(2) 


5.20 


0.300(1) 


945(5) 


0.579(6) 


164(1) 


5.25 


0.284(5) 


945(20) 


0.682(36) 


173(3) 


5.30 


0.260(1) 


973(9) 


0.889(10) 


189(1) 


5.35 


0.233(2) 


959(14) 


1.235(32) 


211(2) 



computational cost constrains us to use only 4^ lattices 
for this study. 

The computational cost increases linearly with the pa- 
rameter N used to define the Fourier transform. For our 
study, we used iV = 12. For each value of (3 we run three 
simulations: n = 0, n = 3, and n = 6. They correspond 
to 0, 1, and 2 baryons in the box. 

Since our volume in lattice units is small, we had 
to use large lattice spacings. We had runs for (3 = 
5.00,5.10,5.15,5.20,5.25,5.30, and 5.35 and we fixed 
K = 0.158. The relevant parameters can be found in 
Table n The lattice spacing and the pion mass are deter- 
mined using standard dynamical action on a 12^ lattice 
for the same values of f3 and k. The lattice spacing was 
determined by using ro scale [23 ■ We note that the pion 
mass varies very little with (3, consequently the quark 
mass is roughly the same in all runs. We also note that 
the quark mass is quite heavy, above the strange quark 
mass. 

For the HMC update, we used the $ algorithm 
made exact by an accept/reject step at the end of each 
trajectory "251. For updating process, we set the length 
of the trajectories to 0.5 with At = 0.01. The HMC 
acceptance rate was very close to 1 since the step length 
was very small. We adjust the number of HMC trajecto- 
ries between two consecutive finite density accept/reject 
steps so that the acceptance rate stays in the range 15% 
to 30%. The relevant information is collected in Table 
im We see that we can get decent acceptance rates even 
when consecutive finite density Metropolis steps are quite 
far apart in configuration space. This allows us to move 
very fast through the configuration space. We collected 
about 100 configurations for each run, separated by 10 
accept/reject steps. 

From an algorithmic point of view, one of the most 
interesting questions is whether or not we have a sign 
problem. To settle this question, we measured the aver- 
age phase a{U) given in Eq. (|27|l . In fact, it is easy to 
prove that the imaginary part of the phase should vanish 
on the ensemble average. It is the real part of this phase 
that carries the signal of a sign problem; if its average is 
close to zero then we have a sign problem. We note here 
that the real part of the phase is ±1, and that the sign 
problem appears when we have an almost equal number 
of configurations of each sign. 
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In Fig. 01 we plot the average of the real part of the 
phase a as a function of the temperature. We note that in 
the deconfined phase, the projected determinant is posi- 
tive most of the time; as we go into the hadronic phase 
the sign starts oscillating. Deep in the hadronic phase, 
the oscillations are more severe at higher density which 
we can see by comparing the case of rt = 6 with n = 3 in 
Fig. ^ However, it is possible that at T < Tc the sign 
average is actually smaller at lower density. This is due 
to the fact that, at this temperature, it is possible to have 
the system in the hadronic phase at low densities and in 
the quark-gluon plasma phase at higher densities. Since 
the oscillations are more severe in the hadronic phase it 
is not surprising that close to and below Tc we would 
have more sign oscillations at lower density. This could 
explain the average sign reversal of n = 3 and n = 6 at 
T = 164 MeV as compared to those at other tempera- 
tures in Fig. ^ 

In Fig. 01 we also see that the sign average drops 
sharply as we go through the transition temperature but 
the rate slows down, as we go deeper in the hadronic 
phase. This slowing down may be due to the fact that 
as we go to lower and lower temperatures, the physical 
volume of the box is also increased and the density de- 
creases. 

In conclusion, it seems that, at least for a 4^ lattice, we 
should be able to investigate the region where T > O.ST^ 
and baryon number ub < 3. We listed here the baryon 
number since we expect that the sign fluctuations are 
going to be determined by this number rather than the 
baryon density. From Table HI we see that the densities 
used in this study are rather large; they range from 2.4 
to 24 times the nuclear matter density. An interesting fu- 
ture direction would be to increase the spatial volume, us- 
ing for example a 6'^ x 4 lattice, while keeping the baryon 
number the same. This would allow us to study densities 
closer to the physically interesting region. If the sign os- 
cillation is really determined by the baryon number not 
volume, we should be able to use this algorithm to ex- 
plore this region. It is also clear that a sign problem will 
appear at baryon numbers larger than the ones employed 
in this study. We also show that the algorithm can be 
efficient in going through the configuration space. For 



TABLE II: Acceptance rates; we list first the number of HMC 
trajectories between two consecutive finite density Metropolis 
steps and then the acceptance rate. 
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FIG. 4: Average phase for n — 3 and n = 6 runs. The n = 
runs have all signs equal to 1 so we did not plot them in this 
figure. 



non-zero density runs, the acceptance rate drops quite 
significantly with the temperature; much smaller number 
of trajectories have to be used between the successive ac- 
cept/reject steps. This may be due to a decrease in the 
autocorrelation time; as we go to smaller values of /3 the 
autocorrelation is expected to decrease. A more detailed 
study is needed to quantify this statement. 



V. PHYSICAL RESULTS 

We turn now toward the physical results. We will 
present measurements of the Polyakov loop, chemical po- 
tential, chiral condensate and the conserved charge. We 
feel compelled to point out that the results presented here 
have large systematic errors. The lattice volume and the 
baryon number are small, consequently the finite size ef- 
fects are going to be important. The lattice spacing is 
very large and the lattice artifacts will be substantial. 
Since we are using Wilson fermions we expect that the 
chiral symmetry is broken quite badly by lattice terms. 
Also, the quark mass is rather heavy. In light of these 
problems, the results presented in this section are inter- 
esting more as proof of concept results. 



A. Polyakov loop 

The most straightforward way to look for a deconfining 
transition is to measure the Polyakov loop. Although the 
average value is expected to vanish due to the sym- 
metry, we can look at the average absolute value. This is 
expected to increase sharply as we go from the confined 
to the deconfined phase. To measure the Polyakov loop 
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FIG. 5: Polyakov loop as a function of temperature. 



we need to fold in the phase 



FIG. 6: A schematic view of the expected QCD phase diagram 
in the temperature-density plane j^. The dotted line, n > 0, 
represents the trajectory in the phase space when we keep the 
baryon number fixed and vary /3. 



1^1) 



(")o 



(36) would produce 



where we denoted with ()q the average over generated en- 
semble. In Fig. [SJ we plot the Polyakov loop for our three 
sets of simulations as a function of temperature. We see 
that a transition occurs somewhere around 170 MeV for 
zero density and, as we increase the density, the transi- 
tion becomes less sharp and moves to lower temperature. 
This picture agrees with the expectations from a study 
with static quarks 's'l, since at large densities the tran- 
sition is expected to be first order and, as a result, the 
system will develop a coexistence region. To visualize 
this, we plot in Fig. 0the expected phase diagram in the 
temperature-density plane. The main difference from the 
picture in the temperature - chemical potential plane (see 
Fig. nj is that the first order transition line is split; we 
have now a line that borders the pure hadronic phase 
and another that borders the pure quark-gluon plasma 
phase. In between them, we have a coexistence region 
characteristic of a first order phase transition. As we go 
through this region, we expect a more pronounced slope 
in \P\. In the infinite volume limit, we expect that the 
slope will change abruptly as we go through the phase 
boundaries but we will not see an abrupt jump in our 
finite density study. 



B. Chemical potential 

In order to compare our results with the results in the 
grand canonical ensemble, we need to measure the chem- 
ical potential. The thermodynamic definition 



1 



Zc{V,T,n) 

N-l 



(38) 



l5]0,e-»*^-detM2(C/^J 



3=0 



where (3 = l/ksT. There are a number of problems with 
this definition. Firstly, the partition function is symmet- 
ric under the transformation 0^ — s- + 2'k. This is not 
true for this definition of the chemical potential and then 
we can ask why would the chemical potential depend 
on our choice of (f>j. Secondly, the chemical potential 
defined above is the quark chemical potential; it mea- 
sures the response of the system when one more quark 
is introduced in the system. If we follow the same logic 
and measure the baryon chemical potential we find that 
A*b("-b) = (*30)3„^ = 3/.t(3ris). Thus, it seems that the 
response to introducing a baryon in the system is linearly 
related to the quark chemical potential. While this might 
be true in the deconfined phase, it is clearly not so in the 
confined phase. The cost of introducing one quark in an 
empty box should be infinite, whereas we expect that the 
cost of introducing a baryon should be finite. To address 
these shortcomings, we "discretize" the derivative and 
define the chemical potential 



^(n) 



F{n + 1) - F{n) 
(n + I) — n 



= F{n + 1) - F{r 



(39) 



We see that defined as above, the chemical potential mea- 
sures the increase in the free energy as we add a quark 
to the system. We find then 



dF{V, T, n) 
dn 



1 d\nZc{V,T, n) 



(37) 
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FIG. 7: Baryon chemical potential as a function of the tem- 
perature. We did not include the data for n = 6 since it 
replicates the results for n = 3. 
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Similarly, for the baryon chemical potential we find 

1 



(41) 



With these new definitions, the quark and baryon chemi- 
cal potentials are no longer linearly related and they also 
satisfy the same symmetries as the partition function. 
Moreover, since the partition function for a system with 
a number of quarks that is not a multiple of 3 vanishes 
when we are in the confined phase, we have 



1 Zc{3n+1) 
u(3n) = — -in — = = +00, 



(42) 



which is exactly what we expect. 

We note that the chemical potential as defined above 
has certain symmetries; since we used Zc for our defini- 
tions we have: 



ju(n) = fi{n + N), 
fJ^Bins) = n{nB + N/3), 



(43) 



the second equality holds when is a multiple of 3. From 
charge conjugation symmetry, we infer that for n, > 



n{-n) = ~fi{n - 1) 
M-b(-?ib) -iJ^sinB - 1) 



(44) 



In Fig. 13 we plot the baryon chemical potential 
as a function of temperature. We see that, as we go 
through the phase transition, the chemical potential 
drops sharply. This is due to the fact that new degrees 



of freedom become available and the entropy of the sys- 
tem increases. We notice that in the confined phase, the 
chemical potential doesn't change much as we increase 
the density, whereas in the deconfined region the chem- 
ical potential is larger as the density increases. These 
findings are consistent with the results of Kratochvila 
and de Forcrand |^ We would also like to point 

out that since in our simulations we used N = 12, we 
can show, using the symmetries of the chemical potential 
listed above, that /U_b(2) = — /xb(1)- This is why we plot 
only the curves for = and 71^ = 1- 

Computing the chemical potential allows us not only 
to connect our results to those from grand canonical sim- 
ulations, but also to determine the shape of the phase 
boundary. To see this, we follow an argument by Kra- 
tochvila and de Forcrand They start by noticing 
that the chemical potential in the hadronic phase seems 
independent of the baryon number. Based on this obser- 
vation they build a simple model where the free energy 
is proportional to the baryon number, F(nB) = ^oI'^^bI- 
The coefficient /xq is the value of the chemical potential 
measured; they show that in this model /zq is just the 
critical chemical potential. Consequently, at any tem- 
perature T if we find /i to be independent of we have 
determined Hc{T). 

The argument above can be generalized. In a physical 
system, the chemical potential is expected to vary from 
small values at small densities to arbitrarly large values as 
the density goes to infinity. However, it can be argued on 
general grounds that when the chemical potential stays 
the same for a range of baryon numbers we are at a phase 
transition. In the thermodynamic limit, the free energy 
is a convex function of the baryon number and thus 



dnr 



> 0. 



(45) 



It is then expected that the chemical potential be an in- 
creasing function of baryon number; it will flatten only 
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as we go through the coexistence region of a first or- 
der phase transition. We see then that we don't really 
need the free energy to be linear in the baryon number; 
when the second derivative vanishes we are at the phase 
boundary. 

In Fig. [71 we see that the chemical potential curves for 
different baryon numbers overlap, at least at the level of 
the error bars, for temperatures lower than 170 MeV. By 
the above argument this part of the curve represents the 
phase boundary and we can use it to plot in Fig. |Slthe 
phase boundary in the (T, /i) plane. 

Before we move on, we would like to point out that, 
although not explicit in the notation we used above, the 
chemical potential, as defined in Eq. (JHl depends also on 
the volume and the temperature of the system. In the 
termodynamic limit the chemical potential, /2, depends 
only on temperature and density: 
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FIG. 9: The quark condensate (^tptpj in lattice units as a 
function of temperature. 



fl{p,T)= \im fi{pV,V,T). (46) 

It is then more relevant to think of our measured chemi- 
cal potential as being defined at a given density. At small 
baryon number there is an ambiguity as to which density 
to assign to a particular measurement since the chemi- 
cal potential is defined to be the difference of the free 
energies at ub + ^ and ns baryon numbers. For large 
values of this does not make much of a difference. 
Most naturally we should think that the measurement is 
performed a.t ns + ^ and treat ^{nB,V,T) as an approx- 
imation for /i(— ^i^,T). For example /i(nB — 0,V,T) 
should be thought as approximating p,{-^,T). It is then 
clear that 24] limy^oo Kns = 0,V,T) = m(0,T) = 0, 
from Eq. 1441 which conforms with expectation. But this 
is not the relevant limit; the limit of interest is that of Eq. 

In other words the density should be kept fixed when 
approaching the thermodynamic limit. Note also that 
using this convention we can show using Eq. 1441 that our 
approximation for the chemical potential becomes sym- 
metric in density, i.e. j2{—p) = —p,{p). 

Another interesting point is that since the chemical po- 
tential p,{nB) should decrease as the volume is increased 
(the density decreases) it would seem that the phase 
boundary constructed using the reasoning we presented 
above will shift. This is not true: the argument rests on 
the fact that the chemical potential stays the same as we 
increase the baryon number; we understand that to be a 
consequence of the fact that we measure the chemical po- 
tential at densities in the phase coexistence region. As we 
increase the volume, the chemical potential will start to 
drop only when we get out of the coexistence region but 
by then it will no longer be independent of hb- To get 
back to a chemical potential that is independent of hb 
we have to increase the baryon number until the density 
again is in the coexistence region. Thus the new phase 
boundary we get at different volumes should be the same 
(up to finite volume corrections). 



C. Quark condensate 

As we cross over from the hadronic phase to deconfined 
phase, we also expect to restore the chiral symmetry. 
There is ample empirical evidence that the deconfining 
phase transition and the chiral symmetry restoration oc- 
cur at almost the same temperatures. As far as we know 
there is no theoretical explanation of this fact, thus it is 
interesting to see whether this remains true at finite den- 
sity. For this purpose, we measure the chiral condensate 
(^ipi/j'). For fermionic observables, we need not only fold 
in the phase a{U); we also need to perform a separate 
Fourier transform. For an arbitrary fermionic bilinear 
^r?/;, where T is some spinor matrix, we have 

X J ViPV^e-^f^^*^''''''''^r^ (47) 

= /|;'feif,_,T„._.,rM-.)\, 

where the factor of 2 comes from using two degenerate 
fiavors and we defined 

Tr„rAf-i = jjY. TrrM([/^J-i, (48) 

the n*'' Fourier component of the trace. Note that when 
computing a fermionic observable, we have contributions 
not only from the O*'* component det„M^TrorM~^, but 
also from the parts of the propagator that wrap around 
the lattice in the time direction. More importantly, de- 
terminant sectors other than det„M^ become relevant. 
To look for the chiral restoration phase transition, we 
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measure the chiral condensate 

(49) 

where Nf = 2 is the number of flavors and A^4 is the 
lattice four volume. In Fig. O we plot the quark con- 
densate as measured in our simulations. We note that 
as we go through the phase transition, the quark con- 
densate gets smaller. The slope changes in the proper 
temperature range and the transition temperature seem 
to decrease with increasing baryon number. Unfortu- 
nately, it is not easy to attach much meaning here, since 
with Wilson fermions the chiral symmetry is broken by 
lattice artifacts and the quark condensate receives large 
contributions from these artifacts. Also, the quark mass 
we employed is very large so the explicit breaking of the 
chiral symmetry is probably large enough to prevent one 
from seeing any signal in the chiral condensate. 



D. Conserved charge 

Finally, we turn our attention to the conserved charge. 
While in the grand canonical ensemble measuring the 
conserved charge, 

Q{t) = ^kYP{x)Ua{x){1 - 74)7^(2; + i) 

X 

-i^{x + i)ul{x){l+-fi)i^{x)], (50) 

helps in measuring the average number of particles in 
the box, there seems to be little point in measuring the 
conserved charge in the canonical ensemble. In fact, we 
can prove that if you are to use the true partition function 
Zc given in Eq. the charge should be equal to the 

number of fermions that we put in the box, configuration 
by configuration. However, since we are simulating an 
approximation of the partition function, Zc, we can use 
the conserved charge to check whether our assumption 
that Zc ~ Zc is true. It is easy to show that 



Conserved charge 



m)) 



Zc(n) 



Y.^{n + mN)Zc{n + mN) 



(51) 



We see then that the deviation from the expected number 
of quarks would quantify how much mixing of different 
quark sectors we have. 

In our simulations, we used = 12 and n = 0, 3 
and 6. For the n — Q and n — & simulations, we 
can prove that the conserved charge will be zero. In 
the n = case, this is due to the charge conjugation 
symmetry Zciji) = Zc{—n). For the n = 6 simula- 
tion, this is due to the fact that for every number of 
the form 6 -I- 12to there is another integer m' such that 
6 -|- 12m' = — (6 -|- 12m); plugging this in the expression 
above we get {Q)n^Q = 0. The only non-trivial case is 
when n = 3, which we plotted in Fig. ^| We see that 
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FIG. 10: Total number of particles in the box for the n — 3 
simulations as a function of temperature. 



as the quarks become deconfined, the chemical poten- 
tial drops and the mixing with the other sectors becomes 
more important. However, even for large temperatures, 
T ~ 200 MeV, the mixing is only about one percent. 
This implies that even the nearest sector, Zc(n = 9), is 
greatly suppressed, i.e. Zc{n = 9)/Zc{n = 3) ~ 0.01. 
The mixing is small mainly because the chemical poten- 
tial is large; the mixing will get worse when we decrease 
the fermion mass. 

In conclusion, in this section we show that we can 
see the expected deconfining transition in the Polyakov 
loop, that the chemical potential drops as we go from 
the hadronic phase to the quark-gluon plasma phase, and 
that there is some hint of the transition even in the chiral 
condensate. Using the conserved charge, we have checked 
that the approximation we made, employing a discrete 
Fourier transform rather than a continuous one, is valid. 



VI. CONCLUSIONS AND OUTLOOK 

In this paper, we show that the canonical partition 
function can be used to investigate the phase structure 
of QCD at finite temperature and non-zero density. The 
algorithm we employed allows us to investigate densities 
much higher that those available by other methods. Sign 
fluctuations limit our ability to reach very low temper- 
atures or much larger densities than the ones explored 
here. We also checked that the discrete Fourier approx- 
imation to the canonical partition function introduces 
only minimal deviations. 

The physical picture that emerges from our simulations 
is consistent with expectations. The Polyakov loop, the 
chemical potential and the quark condensate show signs 
of a transition around T ~ 170 MeV. The quark con- 
densate does not vanish but we need to employ smaller 
masses and finer lattices to reasonably expect a clear sig- 
nal of chiral symmetry restoration. Another route is to 
employ a more sophisticated definition for the chiral con- 
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densate, involving perhaps some form of subtraction, or 
use chiral fermions. 

In the future, we would like to locate points on the 
phase transition line. For this, we need to get to lower 
temperatures and densities. While we might be limited in 
reaching lower temperatures, we should be able to reach 
lower densities; all we need is to move to larger volumes. 
However, since we have to use larger lattices we need to 
use an estimator for the determinant. We should point 
out that the method used to generate the ensemble has 
no bearing on whether we have a sign problem or not, it 
is an intrinsic property of the ensemble. Consequently, 
the sign oscillations stay the same even when we employ 
the determinant estimator. The only thing that is going 
to change is the acceptance rate. We anticipate that this 
should not be a problem since the acceptance rate is very 
good for rather large HMC trajectory lengths. However, 
this need to be studied further. 

Before we conclude, we would like to emphasize that, 
even if it proves that it is not feasible to reach lower 
temperatures, this approach is valuable since it permits 
the study at the phase diagram at temperatures close 
to Tc and rather large densities. We will then be able 
to determine various points on the phase transition line. 
Much effort is put nowadays on determining this line. 



and as we pointed out in the beginning, the methods 
used today need to be checked for reliability. To stress 
this point, we plot in Fig. |H1 next to our phase transi- 
tion line, the second order phase transition point as de- 
termined by Fodor and Katz j3||. Their simulations use 
different quark masses, but the shape of the transition 
line is expected to change very little. Although our er- 
ror bars are rather large, the plot suggests the possibility 
of a discrepancy. This has also been noted in and 
a possible explanation is provided in These results 
seem to indicate a possible overlap problem. It is then 
imperative that new simulations are carried out to check 
the validity of this important result. A future study that 
employs a determinant estimator will allow us to collect 
better statistics and hopefully will settle this question. 
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